Calidad del aire en Madrid, 2018

Objetivo

Hacer una primera aproximación a los datos disponibles sobre la contaminación en el municipio de Madrid.

Vamos a analizar los datos de NO₂ de 2018 siguiendo estos pasos:

• Descarga y extracción de los datos

• Limpieza de de los datos

• Análisis univariante de la serie temporal

• Enriquecimiento de los datos con otras fuentes

• Análisis multivariante de los datos

• Planteamiento de un modelo predictivo

Las librerías que vamos a usar son:

pkgs <- c("lubridate", "data.table", "ggplot2", "prophet", "ggExtra", "leaflet", "magrittr", "dplyr", "leaflet", "gplots", "xts")
lapply(pkgs, library, character.only = TRUE, quietly = T)

Descarga y extracción de los datos

Los datos de calidad del aire se han descargado desde el portal de datos abiertos de Madrid y se han almacendo descomprimidos dentro de la carpeta Anio201810 en data_input.

Vamos a unificar la información de los ficheros de la calidad del aire en un único data frame llamado raw_air_quality.

path_files_air <- "data_input/Anio201810/"
file_list <- list.files(path = path_files_air, pattern = ".csv")

raw_air_quality <- data.frame()
for (file in 1:length(file_list)) {
    path_file <- paste0(path_files_air, file_list[file])
    raw_file <- read.csv(path_file, header = TRUE, sep=";")
    raw_air_quality <- rbind(raw_air_quality, raw_file)
}

rm(raw_file)
dim(raw_air_quality)

## [1] 54831    56

Cada fila hace referencia a todas las medidas tomadas hora a hora durante el mismo día en una misma estación.

names(raw_air_quality)

##  [1] "PROVINCIA"      "MUNICIPIO"      "ESTACION"       "MAGNITUD"      
##  [5] "PUNTO_MUESTREO" "ANO"            "MES"            "DIA"           
##  [9] "H01"            "V01"            "H02"            "V02"           
## [13] "H03"            "V03"            "H04"            "V04"           
## [17] "H05"            "V05"            "H06"            "V06"           
## [21] "H07"            "V07"            "H08"            "V08"           
## [25] "H09"            "V09"            "H10"            "V10"           
## [29] "H11"            "V11"            "H12"            "V12"           
## [33] "H13"            "V13"            "H14"            "V14"           
## [37] "H15"            "V15"            "H16"            "V16"           
## [41] "H17"            "V17"            "H18"            "V18"           
## [45] "H19"            "V19"            "H20"            "V20"           
## [49] "H21"            "V21"            "H22"            "V22"           
## [53] "H23"            "V23"            "H24"            "V24"

Transformación de los datos

Vamos a trabajar sólo con los datos de NO₂, cuya magnitud está codificada como 8.

air_quality <- raw_air_quality[raw_air_quality$MAGNITUD == 8, ]

Para simplificar posteriores análisis, añadimos una nueva columna con el código de las estaciones que corresponde con los 8 primeros dígitos del campo PUNTO_MUESTREO. Es la concatenación del los identificadores de provincia más municipio junto con el de la estación.

station <- as.numeric(substr(air_quality$PUNTO_MUESTREO, 1, 8))

Reajustemos el formato para que cada línea contenga sólo datos verificados, por estación, día y hora.

Vamos a usar un data frame, air_quality_by_days, con la información de estación y fecha que utilizaremos como base para desglosar por fila los valores de NO₂ y validez, por cada hora del día.

air_quality_by_days <- data.frame(
    station, year = air_quality$ANO,
    month = ifelse(air_quality$MES < 10, paste0("0", air_quality$MES), air_quality$MES),
    day = air_quality$DIA,
    stringsAsFactors = FALSE
)
str(air_quality_by_days)

## 'data.frame':    8750 obs. of  4 variables:
##  $ station: num  28079004 28079004 28079004 28079004 28079004 ...
##  $ year   : int  2018 2018 2018 2018 2018 2018 2018 2018 2018 2018 ...
##  $ month  : chr  "04" "04" "04" "04" ...
##  $ day    : int  1 2 3 4 5 6 7 8 9 10 ...

Nota: Resto una unidad al índice de hora para que la hora esté en el rango 0 a 23.

air_quality_by_hour <- data.frame()

for (hour in 1:24) {
    air_quality_by_days$hour <- hour - 1
    air_quality_by_days$no2 <- air_quality[ , (2 * hour + 7)]
    air_quality_by_days$verif <- air_quality[ , (2 * hour + 8)]

    air_quality_by_hour <- rbind(air_quality_by_hour, air_quality_by_days)
}

rm(air_quality_by_days)

Nos quedamos sólo con los datos verificados:

air_quality_by_hour$no2[air_quality_by_hour$verif != "V"] <- NA

Y añadimos una columna que contenga la fecha en formato estándar.

air_quality_by_hour$date <- ymd_hms(paste0(
    air_quality_by_hour$year, "-", air_quality_by_hour$month, "-", air_quality_by_hour$day,
    " ", gsub(pattern = "H", replacement = "", air_quality_by_hour$hour), ":00:00"
))

str(air_quality_by_hour)

## 'data.frame':    210000 obs. of  8 variables:
##  $ station: num  28079004 28079004 28079004 28079004 28079004 ...
##  $ year   : int  2018 2018 2018 2018 2018 2018 2018 2018 2018 2018 ...
##  $ month  : chr  "04" "04" "04" "04" ...
##  $ day    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ hour   : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ no2    : num  21 67 14 8 20 71 32 23 13 7 ...
##  $ verif  : Factor w/ 2 levels "N","V": 2 2 2 2 2 2 2 2 2 2 ...
##  $ date   : POSIXct, format: "2018-04-01" "2018-04-02" ...

Análisis de la serie temporal

Exploración de los datos

Primero obtengamos una visión de la tendencia en 2018:

air_quality_by_day <- air_quality_by_hour
air_quality_by_day$date <- as.Date(air_quality_by_day$date, format = "%Y-%m-%d")

dt_no2 <- data.table(air_quality_by_day)[ , .(median_no2 = median(no2, na.rm = T)), by = .(date)]
names(dt_no2) <- c("ds", "y")

ggplot(data = dt_no2, aes(x = ds, y = y)) + geom_line(color = "#1F77B4") +
    geom_smooth(method = 'loess', color = "#1F77B4", linetype = 2) +
    labs(x = "", y = expression(NO[2]))

Y por días de la semana:

air_quality_by_day$weekday <- weekdays(air_quality_by_day$date)
air_quality_by_day$weekday <- factor(
    air_quality_by_day$weekday,
    levels=c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")
)
air_quality_by_day <- air_quality_by_day[order(air_quality_by_day$weekday), ]

ggplot(air_quality_by_day, aes(x = weekday, y = no2)) + geom_boxplot(fill = "#5fa2ce") +
    labs(x = "", y = expression(NO[2]))

## Warning: Removed 919 rows containing non-finite values (stat_boxplot).

Obsevando la componente semanal y la anual podemos llegar a la conclusión de que el nivel de NO₂ baja en días festivos, pero también en los no lectivos.

Como boceto visual de cómo serían los mapas de calor muestro éste de horas frente días del mes (nota: habría que eliminar los datos no verificados y cambiar la leyenda de los ejes para que fuera más fácil de leer).

air_quality <- select(raw_air_quality, -contains("V"))
air_quality_by_day_by_hour <- select(air_quality, 7:ncol(air_quality))
aggdata <- aggregate(
  air_quality_by_day_by_hour,
  by = list(air_quality_by_day_by_hour$DIA),
  FUN = mean, na.rm = TRUE
)
aggdata <- data.matrix(aggdata[ -c(1, 2)])

colheatmap <- c("#1c5998", "#1c73b1", "#3a87b7", "#67add4","#7bc8e2","#cacaca", "#fdab67", "#fd8938","#f06511", "#d74401", "#a33202")
heatmap.2(
  aggdata,
  dendrogram = "none",
  Colv = NA,
  Rowv = NA,
  density.info = "none",
  trace = "none",
  col = colheatmap
)

Como era previsible, en las horas nocturnas, con poco tráfico rodado, también baja el nivel de NO₂. Quedaría pendiente analizar, desglosándolo por meses, por qué a veces se produce una mayor concentración a última hora del día, sobre todo a final de mes.

¿Y si vemos cómo se comportan los valores máximos y mínimos frente a la mediana?

dt_no2 <- data.table(air_quality_by_day)[ , .(
    median_no2 = median(no2, na.rm = T),
    min_no2 = min(no2, na.rm = T),
    max_no2 = max(no2, na.rm = T) ), by = .(date)]

ggplot(dt_no2, aes(date)) +
    geom_line(aes(y = min_no2), color = "#5f9ed1") +
    geom_line(aes(y = max_no2), color = "#fbb04e") +
    geom_line(aes(y = median_no2), color = "#006ba4") +
    labs(x = "", y = expression(NO[2]))

Tanta diferencia entre los valores extremos dentro del mismo día pueden ser debidos a que las mediciones se han tomado en distintas zonas geográficas. Debemos introducir más fuentes de datos para poder seguir investigando.

Forecasting

He probado Prophet, pero en este caso no lo recomiendo pues hace una previsión de NO₂ muy alta para enero de 2019, lo que desvirtúa el resto del análisis. Posiblemente se corregiría si incluyéramos más años en el análisis, o bien podríamos hacer un análisis ARIMA.

dt_no2 <- data.table(air_quality_by_day)[ , .(median_no2 = median(no2, na.rm = T)), by = .(date)]
names(dt_no2) <- c("ds", "y")

m <- prophet(dt_no2, daily.seasonality = TRUE, yearly.seasonality = TRUE )

future <- make_future_dataframe(m, periods = 360)
forecast <- predict(m, future)
plot(m, forecast)

prophet_plot_components(m, forecast)

Otros análisis mulltivariantes que podrían ser interesantes:

• Análisis de ANOVA determinar si existen diferencias estadísticamente significativas de NO₂ entre las horas del día, entre días de la semana, entre días festivos, entre días lectivos, entre días del mes o entre los meses.
• Mapa de calor entre día de la semana y hora.
• Mapa de calor entre día de la semana y día del mes.
• Estudiar los puntos outliers que se ven en el boxplot NO₂ de por semana.

Enriquecimiento de los datos con otras fuentes

Vamos a enriquecer estos datos con la geolocalización de las estaciones de medición, así como la medida de la temperatura.

Utilizaremos la información de las estaciones contenida en el fichero madrid_air_quality_stations.csv de la carpeta data_input, extraído de aquí:

stations <- read.csv("data_input/madrid_air_quality_stations.csv", header = TRUE, sep = ",")
str(stations)

## 'data.frame':    24 obs. of  8 variables:
##  $ station  : int  28079004 28079008 28079011 28079016 28079017 28079018 28079024 28079027 28079035 28079036 ...
##  $ area     : Factor w/ 5 levels "centro","noreste",..: 1 1 1 2 5 5 3 2 1 4 ...
##  $ name     : Factor w/ 24 levels "Arturo Soria",..: 17 10 2 1 24 11 5 3 18 14 ...
##  $ address  : Factor w/ 24 levels "Avd. Betanzos esq. C/ Monforte de Lemos",..: 21 17 6 8 13 15 16 14 22 2 ...
##  $ altitude : int  635 670 708 693 604 630 642 621 659 685 ...
##  $ type     : Factor w/ 3 levels "S","UF","UT": 3 3 3 2 2 2 1 2 2 3 ...
##  $ longitude: num  -3.71 -3.68 -3.68 -3.64 -3.71 ...
##  $ latitude : num  40.4 40.4 40.5 40.4 40.3 ...

df_air <- merge(air_quality_by_hour, stations)
str(df_air)

## 'data.frame':    210000 obs. of  15 variables:
##  $ station  : num  28079004 28079004 28079004 28079004 28079004 ...
##  $ year     : int  2018 2018 2018 2018 2018 2018 2018 2018 2018 2018 ...
##  $ month    : chr  "04" "04" "04" "04" ...
##  $ day      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ hour     : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ no2      : num  21 67 14 8 20 71 32 23 13 7 ...
##  $ verif    : Factor w/ 2 levels "N","V": 2 2 2 2 2 2 2 2 2 2 ...
##  $ date     : POSIXct, format: "2018-04-01 00:00:00" "2018-04-02 00:00:00" ...
##  $ area     : Factor w/ 5 levels "centro","noreste",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ name     : Factor w/ 24 levels "Arturo Soria",..: 17 17 17 17 17 17 17 17 17 17 ...
##  $ address  : Factor w/ 24 levels "Avd. Betanzos esq. C/ Monforte de Lemos",..: 21 21 21 21 21 21 21 21 21 21 ...
##  $ altitude : int  635 635 635 635 635 635 635 635 635 635 ...
##  $ type     : Factor w/ 3 levels "S","UF","UT": 3 3 3 3 3 3 3 3 3 3 ...
##  $ longitude: num  -3.71 -3.71 -3.71 -3.71 -3.71 ...
##  $ latitude : num  40.4 40.4 40.4 40.4 40.4 ...

Utilizaremos la información de temperaturas medias por días contenidas en el archivo madrid_hourly_temperatures_2018.csv de la carpeta data_input, extraído de aquí.

Llamaremos date_day a la columna con la información del día y date a una nueva columna que incluya también la hora:

temperatures <- read.csv("data_input/madrid_hourly_temperatures_2018.csv", header = TRUE, sep = ",")
temperatures$date_day <- ymd(temperatures$date)
temperatures$date <- ymd_hms(paste(temperatures$date_day, temperatures$hour), truncated = 3)
str(temperatures)

## 'data.frame':    8760 obs. of  4 variables:
##  $ date    : POSIXct, format: "2018-01-01 00:00:00" "2018-01-01 01:00:00" ...
##  $ hour    : int  0 1 2 3 4 5 6 7 8 9 ...
##  $ temp    : num  7.3 7.4 6.8 7.1 5.3 3.8 2.9 2.9 2.3 4.2 ...
##  $ date_day: Date, format: "2018-01-01" "2018-01-01" ...

Guardamos toda la información en un único data frame y lo exportamos a un fichero para futuros análisis:

df_air <- merge(df_air, temperatures)
write.csv(df_air, file = "data_output/df_air.csv", row.names = FALSE)
str(df_air)

## 'data.frame':    210000 obs. of  17 variables:
##  $ hour     : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ date     : POSIXct, format: "2018-01-01" "2018-01-01" ...
##  $ station  : num  28079060 28079038 28079008 28079057 28079055 ...
##  $ year     : int  2018 2018 2018 2018 2018 2018 2018 2018 2018 2018 ...
##  $ month    : chr  "01" "01" "01" "01" ...
##  $ day      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ no2      : num  8 13 33 5 9 12 5 32 37 8 ...
##  $ verif    : Factor w/ 2 levels "N","V": 2 2 2 2 2 2 2 2 2 2 ...
##  $ area     : Factor w/ 5 levels "centro","noreste",..: 2 1 1 2 2 1 2 5 5 4 ...
##  $ name     : Factor w/ 24 levels "Arturo Soria",..: 21 7 10 20 22 4 12 11 19 9 ...
##  $ address  : Factor w/ 24 levels "Avd. Betanzos esq. C/ Monforte de Lemos",..: 23 5 17 11 12 1 18 15 24 3 ...
##  $ altitude : int  715 698 670 700 618 674 660 630 604 627 ...
##  $ type     : Factor w/ 3 levels "S","UF","UT": 2 3 3 2 2 3 1 2 3 2 ...
##  $ longitude: num  -3.69 -3.71 -3.68 -3.66 -3.58 ...
##  $ latitude : num  40.5 40.4 40.4 40.5 40.5 ...
##  $ temp     : num  7.3 7.3 7.3 7.3 7.3 7.3 7.3 7.3 7.3 7.3 ...
##  $ date_day : Date, format: "2018-01-01" "2018-01-01" ...

Análisis multivariante de los datos

Al analizar la relación entre el nivel de NO₂ y la temperatura no se observa que exista correlación:

p <- ggplot(df_air, aes(x = temp, y = no2, color = -temp)) +
    geom_point() +
    theme(legend.position = "none") +
    labs(x = "Temperatura", y = expression(NO[2]))
ggMarginal(p, type = "histogram", fill = "#1F77B4")

## Warning: Removed 919 rows containing missing values (geom_point).

Ni aunque lo desglosemos por cada una de las estaciones:

ggplot(df_air, aes(x = temp, y = no2, color = -temp)) +
    geom_point() +
    theme(legend.position = "none") +
    facet_wrap( ~ name) +
    labs(x = "Temperatura", y = expression(NO[2]))

## Warning: Removed 919 rows containing missing values (geom_point).

Lo que sí nos da una pista es que existe mucha variación entre el nivel de NO₂ de las distintas estaciones.

Vamos a comparar los valores medianos haciendo tres grupos según al cuartil en el que estén los datos:

station_media <- df_air %>%
    group_by(name) %>%
    summarise( no2 = median(no2, na.rm = TRUE)) %>%
    arrange(desc(no2))

station_media <- as.data.frame(station_media)
station_media$name <- factor(station_media$name,
                             levels = station_media$name)
station_media <- station_media[order(station_media$name), ]

station_media$level_no2 <- "high"
station_media$level_no2[station_media$no2 < quantile(station_media$no2)[4]] <- "medium"
station_media$level_no2[station_media$no2 < quantile(station_media$no2)[2]] <- "low"

ggplot(station_media, aes(x=name, y = no2, group = factor(level_no2), fill = factor(level_no2))) +
    scale_fill_manual(values = c("#fc7d0b", "#1170aa", "#5fa2ce")) +
    geom_bar(stat = "identity") +
    coord_flip() +
    theme_minimal() +
    theme(legend.position = "none") +
    labs(x = "", y = expression(NO[2]))

Y la relación entre las estaciones y su localización geográfica:

station_media <- merge(station_media, stations)

station_media$color_no2 <- "red"
station_media$color_no2[station_media$no2 < quantile(station_media$no2)[4] ] <- "orange"
station_media$color_no2[station_media$no2 < quantile(station_media$no2)[2] ] <- "green"

icons <- awesomeIcons(
    icon = 'ios-close',
    iconColor = station_media$color_no2,
    library = 'ion',
    markerColor = station_media$color_no2
)

leaflet(station_media) %>% addTiles() %>%
    addAwesomeMarkers(lng = station_media$longitude,
                      lat = station_media$latitude,
                      icon = icons,
                      label = ~paste(name, "-", as.character(no2)))

Podemos apreciar cómo las zonas más céntricas contienen mayor nivel de NO₂. La excepción es la estación que se encuentra situada en el Parque de «El Retiro», que tiene un nivel de 22, menos de la mitad que «Escuelas Aguirre» con 51.

Nota: A la vista de los resultados deberíamos incluir «Castellana» dentro del grupo high.

Planteamiento de modelo de predicción

Variable a predecir: NO₂

Variables independientes:

• Estación se encuentra en zona verde
• Día de la semana
• Día festivo laboral
• Día festivo escolar
• Temperatura
• Velocidad del viento
• Humedad
• Precipitación

Sería deseable tener variables meteorológicas por cada estación por día e incluso por hora.

Además, tener histórico de más años, nos va a ayudaría a determinar la estacionalidad de los datos.

Primero dividiríamos el conjunto de datos en dos grupos: entrenamiento y testeo.

Para evitar sobreajustes podemos personalizar los parámetros del conjunto de entrenamiento con la función trainControl de la librería caret, por ejemplo, haciendo k-fold cross-validation que se repita r veces

library(caret)
control <- trainControl(method = "repeatedcv", number = k, repeats = r)

Podríamos probar con los algoritmos: randomForest, xgboost, glmnet para ver cuál hace una mejor predicción:

library(randomForest)
library(xgboost)
library(glmnet)

fit.rf <- train(no2~. , data, method = "rf", trControl = control, na.action = na.exclude)
fit.xg <- train(no2~. , data, method = "xgbLinear", trControl = control, na.action = na.exclude)
fit.gl <- train(no2~. , data, method = "glmnet", trControl = control, na.action = na.exclude)

results <- resamples(list(RF = fit.rf, XG = fit.xg, GL = fit.gl))

Y, por último, tendríamos que analizar qué modelo nos conviene tomar:

summary(results)
scales <- list(x = list(relation = "free"), y = list(relation = "free"))
bwplot(results, scales = scales)

Aquí hay más ejemplos del análisis para decidir entre los distintos algoritmos.

Resumen

1. Descarga y extracción de los datos.
   • Carga de los datos de calidad del aire de Madrid.
   • Pendiente: Tomar al menos otro año para mejorar los análisis posteriores.
   • Conclusiones: Hay que tratar los datos pues no viene en formato estándar.
2. Limpieza de de los datos.
   • Filtrado el conjunto de datos solo para la magnitud NO₂.
   • Transformación del formato para que cada línea contenga datos verificados por estación, día y hora.
   • Eliminación de los datos no validados.
   • Agregación de una columna con formato de fecha.
3. Análisis univariante de la serie temporal.
   • Exploración visual de los datos de forma anual, semanal, y horas frente días del mes.
   • Análisis predictivo con Prophet.
   • Pendiente:
     • Análisis de ANOVA para ver las diferencias significativas de NO₂ entre las horas del día, entre días de la semana, entre días del mes o entre los meses.
     • Mapa de calor entre día de la semana y hora.
     • Mapa de calor entre día de la semana y día del mes.
     • Estudiar los puntos outliers que se ven en el boxplot NO₂ por semana.
     • Análisis predictivo ARIMA.
   • Conclusiones: El nivel de NO₂ baja en días festivos, en los no lectivos y en las horas nocturnas.
4. Enriquecimiento de los datos con otras fuentes.
   • Tratamiento de los datos para incluir datos de temperaturas y de estaciones.
   • Pendiente: Añadir más fuentes como días festivos y días lectivos.
5. Análisis multivariante de los datos.
   • Análisis de la relación entre nivel de NO₂ con la temperatura y con las distintas estaciones de medición.
   • Pendiente:
     • Mapa de calor animado sobre el mapa de Madrid para ver cómo evoluciona el nivel de NO₂ con el temperatura según las estaciones.
     • Análisis de ANOVA para ver las diferencias significativas de NO₂ entre entre días festivos y otro entre días lectivos.
   • Conclusiones:
     • No se ve relación entre el nivel de NO₂ y la temperatura.
     • Podemos apreciar cómo las zonas más céntricas contienen mayor nivel de NO₂. La excepción es la estación que se encuentra situada en el Parque de «El Retiro», que tiene un nivel de 22, menos de la mitad que «Escuelas Aguirre» con 51.
6. Planteamiento de modelo predictivo.
   • Planteamiento de cómo hacer el modelo.
   • Pendiente: Desarrollarlo.